<!DOCTYPE html>

<html xmlns="http://www.w3.org/1999/xhtml">

<head>

<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />

<meta name="viewport" content="width=device-width, initial-scale=1">

<meta name="author" content="Jared Huling" />

<meta name="date" content="2018-03-25" />

<title>Usage of the oem Package</title>



<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
a.sourceLine { display: inline-block; line-height: 1.25; }
a.sourceLine { pointer-events: none; color: inherit; text-decoration: inherit; }
a.sourceLine:empty { height: 1.2em; position: absolute; }
.sourceCode { overflow: visible; }
code.sourceCode { white-space: pre; position: relative; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
code.sourceCode { white-space: pre-wrap; }
a.sourceLine { text-indent: -1em; padding-left: 1em; }
}
pre.numberSource a.sourceLine
  { position: relative; }
pre.numberSource a.sourceLine:empty
  { position: absolute; }
pre.numberSource a.sourceLine::before
  { content: attr(data-line-number);
    position: absolute; left: -5em; text-align: right; vertical-align: baseline;
    border: none; pointer-events: all;
    -webkit-touch-callout: none; -webkit-user-select: none;
    -khtml-user-select: none; -moz-user-select: none;
    -ms-user-select: none; user-select: none;
    padding: 0 4px; width: 4em;
    color: #aaaaaa;
  }
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa;  padding-left: 4px; }
div.sourceCode
  {  }
@media screen {
a.sourceLine::before { text-decoration: underline; }
}
code span.al { color: #ff0000; font-weight: bold; } /* Alert */
code span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #7d9029; } /* Attribute */
code span.bn { color: #40a070; } /* BaseN */
code span.bu { } /* BuiltIn */
code span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4070a0; } /* Char */
code span.cn { color: #880000; } /* Constant */
code span.co { color: #60a0b0; font-style: italic; } /* Comment */
code span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #ba2121; font-style: italic; } /* Documentation */
code span.dt { color: #902000; } /* DataType */
code span.dv { color: #40a070; } /* DecVal */
code span.er { color: #ff0000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #40a070; } /* Float */
code span.fu { color: #06287e; } /* Function */
code span.im { } /* Import */
code span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #007020; font-weight: bold; } /* Keyword */
code span.op { color: #666666; } /* Operator */
code span.ot { color: #007020; } /* Other */
code span.pp { color: #bc7a00; } /* Preprocessor */
code span.sc { color: #4070a0; } /* SpecialChar */
code span.ss { color: #bb6688; } /* SpecialString */
code span.st { color: #4070a0; } /* String */
code span.va { color: #19177c; } /* Variable */
code span.vs { color: #4070a0; } /* VerbatimString */
code span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
</style>



<style type="text/css">

body {
background-color: #fff;
margin: 1em auto;
max-width: 800px;
overflow: visible;
padding-left: 2em;
padding-right: 2em;
font-family: "Palatino Linotype", "Book Antiqua", Palatino, serif;
font-size: 16px;
line-height: 1.35;
}
#header {
text-align: center;
}
#TOC {
clear: both;
margin: 0 0 10px 10px;
padding: 4px;
width: 780px;
border: 1px solid #CCCCCC;
border-radius: 5px;
background-color: #f6f6f6;
font-size: 14px;
line-height: 1.3;
}
#TOC .toctitle {
font-weight: bold;
font-size: 16px;
margin-left: 5px;
}
#TOC ul {
padding-left: 40px;
margin-left: -1.5em;
margin-top: 5px;
margin-bottom: 5px;
}
#TOC ul ul {
margin-left: -2em;
}
#TOC li {
line-height: 16px;
}
table {
margin: 1em auto;
border-width: 1px;
border-color: #DDDDDD;
border-style: outset;
border-collapse: collapse;
}
table th {
border-width: 2px;
padding: 5px;
border-style: inset;
}
table td {
border-width: 1px;
border-style: inset;
line-height: 18px;
padding: 5px 5px;
}
table, table th, table td {
border-left-style: none;
border-right-style: none;
}
table thead, table tr.even {
background-color: #f7f7f7;
}
p {
margin: 0.5em 0;
}
blockquote {
background-color: #f6f6f6;
padding: 0.25em 0.75em;
}
hr {
border-style: solid;
border: none;
border-top: 1px solid #777;
margin: 28px 0;
}
dl {
margin-left: 0;
}
dl dd {
margin-bottom: 13px;
margin-left: 13px;
}
dl dt {
font-weight: bold;
}
ul {
margin-top: 0;
}
ul li {
list-style: circle outside;
}
ul ul {
margin-bottom: 0;
}
pre, code {
background-color: #f7f7f7;
border-radius: 5px;
color: #333;
white-space: pre-wrap; 
}
pre {
border-radius: 5px;
margin: 5px 0px 10px 0px;
padding: 10px;
}
pre:not([class]) {
background-color: #f7f7f7;
}
code {
font-family: "Lucida Console", Monaco, Consolas, 'Courier New', monospace;
font-size: 85%;
}
p > code, li > code {
padding: 2px 0px;
}
div.figure {
text-align: center;
}
img {
background-color: #FFFFFF;
padding: 2px;
border: 1px solid #DDDDDD;
border-radius: 3px;
border: 1px solid #CCCCCC;
margin: 0 5px;
}
h1 {
margin-top: 25px;
font-size: 35px;
line-height: 40px;
}
h2 {
border-bottom: 2px solid #f7f7f7;
padding-top: 10px;
padding-bottom: 2px;
font-size: 145%;
}
h3 {
border-bottom: 1px solid #f7f7f7;
padding-top: 10px;
font-size: 120%;
}
h4 {
margin-left: 8px;
font-size: 105%;
}
em {
font-style: normal;
}
h5, h6 {
border-bottom: 1px solid #ccc;
font-size: 105%;
}
a {
color: #123d79;
text-decoration: none;
}
a:hover {
color: #007CC3; }
a:visited {
color: #581858; }
a:visited:hover {
color: #007CC3; }

code > span.kw { color: #4271ae; } 
code > span.dt { color: #c82829; } 
code > span.dv { color: #f5871f; } 
code > span.bn { color: #718c00; } 
code > span.fl { color: #718c00; } 
code > span.ch { color: #718c00; } 
code > span.st { color: #718c00; } 
code > span.co { color: #8e908c; } 
code > span.ot { color: #4d4d4c; } 
code > span.al { color: #ff0000; } 
code > span.fu { color: #4271ae; } 
code > span.er { color: #a61717; } 

img {
display: block;
margin: 0 auto;
}
</style>

</head>

<body>




<h1 class="title toc-ignore">Usage of the oem Package</h1>
<h4 class="author"><em>Jared Huling</em></h4>
<h4 class="date"><em>2018-03-25</em></h4>


<div id="TOC">
<ul>
<li><a href="#introduction"><span class="toc-section-number">1</span> Introduction</a></li>
<li><a href="#installation"><span class="toc-section-number">2</span> Installation</a></li>
<li><a href="#quick-start"><span class="toc-section-number">3</span> Quick Start</a></li>
<li><a href="#key-features"><span class="toc-section-number">4</span> Key Features</a><ul>
<li><a href="#available-functions"><span class="toc-section-number">4.1</span> Available functions</a></li>
<li><a href="#available-penalties"><span class="toc-section-number">4.2</span> Available Penalties</a></li>
<li><a href="#available-model-families"><span class="toc-section-number">4.3</span> Available Model Families</a></li>
</ul></li>
<li><a href="#fitting-multiple-penalties-at-once"><span class="toc-section-number">5</span> Fitting multiple penalties at once</a><ul>
<li><a href="#timing-comparison"><span class="toc-section-number">5.1</span> Timing Comparison</a><ul>
<li><a href="#linear-regression"><span class="toc-section-number">5.1.1</span> Linear Regression</a></li>
<li><a href="#logistic-regression"><span class="toc-section-number">5.1.2</span> Logistic Regression</a></li>
</ul></li>
</ul></li>
<li><a href="#cross-validation"><span class="toc-section-number">6</span> Cross Validation</a><ul>
<li><a href="#extremely-fast-cross-validation-for-linear-models"><span class="toc-section-number">6.1</span> Extremely Fast Cross Validation for Linear Models</a></li>
<li><a href="#evaluation-metrics"><span class="toc-section-number">6.2</span> Evaluation Metrics</a><ul>
<li><a href="#misclassification-rate"><span class="toc-section-number">6.2.1</span> Misclassification Rate</a></li>
<li><a href="#area-under-the-roc-curve"><span class="toc-section-number">6.2.2</span> Area Under the ROC Curve</a></li>
</ul></li>
</ul></li>
<li><a href="#methods-for-very-large-scale-problems"><span class="toc-section-number">7</span> Methods for Very Large Scale Problems</a><ul>
<li><a href="#oem-with-precomputed-xtx-xty-for-linear-models"><span class="toc-section-number">7.1</span> OEM with Precomputed <span class="math inline">\(X^TX\)</span>, <span class="math inline">\(X^TY\)</span> for Linear Models</a></li>
<li><a href="#out-of-memory-computation"><span class="toc-section-number">7.2</span> Out-of-memory Computation</a></li>
</ul></li>
<li><a href="#other-features"><span class="toc-section-number">8</span> Other Features</a><ul>
<li><a href="#parallelization-via-openmp"><span class="toc-section-number">8.1</span> Parallelization via OpenMP</a></li>
<li><a href="#penalty-adjustment"><span class="toc-section-number">8.2</span> Penalty Adjustment</a></li>
</ul></li>
<li><a href="#more-information"><span class="toc-section-number">9</span> More Information</a></li>
<li><a href="#references">References</a></li>
</ul>
</div>

<div id="introduction" class="section level1">
<h1><span class="header-section-number">1</span> Introduction</h1>
<p><code>oem</code> is a package for the estimation of various penalized regression models using the oem algorithm of <span class="citation">(Xiong et al. 2016)</span>. The focus of <code>oem</code> is to provide high performance computation for big <strong>tall</strong> data. Many applications not only have a large number of variables, but a vast number of observations; <code>oem</code> is designed to perform well in these settings.</p>
<ul>
<li>Fast computation for big <strong>tall</strong> data</li>
<li>Efficient computation for computation for multiple penalties simultaneously</li>
<li>Efficient cross-validation</li>
</ul>
</div>
<div id="installation" class="section level1">
<h1><span class="header-section-number">2</span> Installation</h1>
<p>The simplest way to install <code>oem</code> is via the CRAN repositories as the following:</p>
<p>To install the development version, first install the <code>devtools</code> package</p>
<div class="sourceCode" id="cb1"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb1-1" data-line-number="1"><span class="kw">install.packages</span>(<span class="st">&quot;devtools&quot;</span>, <span class="dt">repos =</span> <span class="st">&quot;http://cran.us.r-project.org&quot;</span>)</a></code></pre></div>
<p>and then installl <code>oem</code> via the <code>install_github</code> function</p>
<div class="sourceCode" id="cb2"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb2-1" data-line-number="1"><span class="kw">library</span>(devtools)</a>
<a class="sourceLine" id="cb2-2" data-line-number="2"><span class="kw">install_github</span>(<span class="st">&quot;jaredhuling/oem&quot;</span>)</a></code></pre></div>
</div>
<div id="quick-start" class="section level1">
<h1><span class="header-section-number">3</span> Quick Start</h1>
<p>First load <code>oem</code></p>
<div class="sourceCode" id="cb3"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb3-1" data-line-number="1"><span class="kw">library</span>(oem)</a></code></pre></div>
<p>Simulate data</p>
<div class="sourceCode" id="cb4"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb4-1" data-line-number="1">nobs  &lt;-<span class="st"> </span><span class="fl">1e4</span></a>
<a class="sourceLine" id="cb4-2" data-line-number="2">nvars &lt;-<span class="st"> </span><span class="dv">100</span></a>
<a class="sourceLine" id="cb4-3" data-line-number="3">x &lt;-<span class="st"> </span><span class="kw">matrix</span>(<span class="kw">rnorm</span>(nobs <span class="op">*</span><span class="st"> </span>nvars), <span class="dt">ncol =</span> nvars)</a>
<a class="sourceLine" id="cb4-4" data-line-number="4">y &lt;-<span class="st"> </span><span class="kw">drop</span>(x <span class="op">%*%</span><span class="st"> </span><span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">-0.5</span>, <span class="fl">-0.5</span>, <span class="dv">1</span>, <span class="kw">rep</span>(<span class="dv">0</span>, nvars <span class="op">-</span><span class="st"> </span><span class="dv">5</span>))) <span class="op">+</span><span class="st"> </span><span class="kw">rnorm</span>(nobs, <span class="dt">sd =</span> <span class="dv">4</span>)</a></code></pre></div>
<p>Fit a penalized regression model using the <code>oem</code> function</p>
<div class="sourceCode" id="cb5"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb5-1" data-line-number="1">fit1 &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> x, <span class="dt">y =</span> y, <span class="dt">penalty =</span> <span class="st">&quot;lasso&quot;</span>)</a></code></pre></div>
<p>Plot the solution path</p>
<div class="sourceCode" id="cb6"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb6-1" data-line-number="1"><span class="kw">plot</span>(fit1)</a></code></pre></div>
<p><img src="" /><!-- --></p>
</div>
<div id="key-features" class="section level1">
<h1><span class="header-section-number">4</span> Key Features</h1>
<div id="available-functions" class="section level2">
<h2><span class="header-section-number">4.1</span> Available functions</h2>
<table>
<thead>
<tr class="header">
<th>Function Name</th>
<th>Functionality</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td><code>oem()</code></td>
<td>Main fitting function</td>
</tr>
<tr class="even">
<td><code>oem.xtx()</code></td>
<td>Fitting function for precomputed <span class="math inline">\(X'X,X'y\)</span></td>
</tr>
<tr class="odd">
<td><code>big.oem()</code></td>
<td>Fitting function for <code>big.matrix()</code> objects</td>
</tr>
<tr class="even">
<td><code>summary.oemfit()</code></td>
<td>Summary for <code>oem</code> objects</td>
</tr>
<tr class="odd">
<td><code>predict.oemfit()</code></td>
<td>Prediction for <code>oem</code> objects</td>
</tr>
<tr class="even">
<td><code>plot.oemfit()</code></td>
<td>Plotting for <code>oem</code> objects</td>
</tr>
<tr class="odd">
<td><code>logLik.oemfit()</code></td>
<td>log Likelihood for <code>oem</code> objects</td>
</tr>
<tr class="even">
<td><code>cv.oem()</code></td>
<td>Cross-validation function</td>
</tr>
<tr class="odd">
<td><code>xval.oem()</code></td>
<td>Fast cross-validation for linear models</td>
</tr>
<tr class="even">
<td><code>summary.cv.oem()</code></td>
<td>Summary for <code>cv.oem</code> objects</td>
</tr>
<tr class="odd">
<td><code>predict.cv.oem()</code></td>
<td>Prediction for <code>cv.oem</code> objects</td>
</tr>
<tr class="even">
<td><code>plot.cv.oem()</code></td>
<td>Plotting for <code>cv.oem</code> objects</td>
</tr>
<tr class="odd">
<td><code>logLik.cv.oem()</code></td>
<td>log Likelihood for cv.oem objects</td>
</tr>
</tbody>
</table>
</div>
<div id="available-penalties" class="section level2">
<h2><span class="header-section-number">4.2</span> Available Penalties</h2>
<table>
<colgroup>
<col width="30%"></col>
<col width="30%"></col>
<col width="40%"></col>
</colgroup>
<thead>
<tr class="header">
<th>Penalty</th>
<th>Option Name</th>
<th>Penalty Form</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>Lasso</td>
<td><code>lasso</code></td>
<td><span class="math inline">\(\lambda \sum_{j = 1}^pd_j|\beta_j|\)</span></td>
</tr>
<tr class="even">
<td>Elastic Net</td>
<td><code>elastic.net</code></td>
<td><span class="math inline">\(\lambda \sum_{j = 1}^p\alpha d_j|\beta_j| + \frac{1}{2}(1 - \alpha)\lambda \sum_{j = 1}^pd_j\beta_j^2\)</span></td>
</tr>
<tr class="odd">
<td>MCP</td>
<td><code>mcp</code></td>
<td><span class="math inline">\(\lambda \sum_{j = 1}^pd_j \int_0^{\beta_j}(1 - x/(\gamma\lambda d_j))_+\mathrm{d}x\)</span></td>
</tr>
<tr class="even">
<td>SCAD</td>
<td><code>scad</code></td>
<td><span class="math inline">\(\sum_{j = 1}^p p^{SCAD}_{\lambda d_j,\gamma}(\beta_j)\)</span></td>
</tr>
<tr class="odd">
<td>Group Lasso</td>
<td><code>grp.lasso</code></td>
<td><span class="math inline">\(\lambda \sum_{k = 1}^Gd_k\sqrt{\sum_{j \in g_k}\beta_j^2}\)</span></td>
</tr>
<tr class="even">
<td>Group MCP</td>
<td><code>grp.mcp</code></td>
<td><span class="math inline">\(\lambda \sum_{k = 1}^Gp^{MCP}_{\lambda d_k,\gamma}(||\boldsymbol\beta_{g_k}||_2)\)</span></td>
</tr>
<tr class="odd">
<td>Group SCAD</td>
<td><code>grp.scad</code></td>
<td><span class="math inline">\(\lambda \sum_{k = 1}^Gp^{SCAD}_{\lambda d_k,\gamma}(||\boldsymbol\beta_{g_k}||_2)\)</span></td>
</tr>
<tr class="even">
<td>Sparse Group Lasso</td>
<td><code>sparse.grp.lasso</code></td>
<td><span class="math inline">\(\lambda \alpha\sum_{j = 1}^pd_j|\beta_j| + \lambda (1-\alpha)\sum_{k = 1}^Gd_k\sqrt{\sum_{j \in g_k}\beta_j^2}\)</span></td>
</tr>
</tbody>
</table>
<p>where <span class="math inline">\(||\boldsymbol\beta_{g_k}||_2 = \sqrt{\sum_{j \in g_k}\beta_j^2}\)</span>, <span class="math display">\[ p_{\lambda, \gamma}^{SCAD}(\beta) =   \left\{
\begin{array}{ll}
      \lambda|\beta| &amp; |\beta| \leq \lambda ; \\
      -\frac{|\beta|^2 - 2\gamma\lambda|\beta| + \lambda^2}{2(\gamma - 1)} &amp; \lambda &lt; |\beta| \leq \gamma\lambda ; \\
      \frac{(\gamma + 1)\lambda^2}{2} &amp; |\beta| &gt; \gamma\lambda, \\
\end{array} 
\right. \]</span></p>
<p>and</p>
<p><span class="math display">\[ p_{\lambda, \gamma}^{MCP}(\beta) =  \lambda \int_0^{\beta}(1 - x/(\gamma\lambda ))_+\mathrm{d}x = \\ \left\{
\begin{array}{ll}
      -\lambda (|\beta| - \frac{\beta^2} {2 \lambda\gamma}) &amp; |\beta| \leq \gamma\lambda ; \\
      \frac{ \lambda^2\gamma}{2} &amp; |\beta| &gt; \gamma\lambda. \\
\end{array} 
\right. \]</span></p>
<p>Any penalty with <code>.net</code> at the end of its name has a ridge term of <span class="math inline">\(\frac{1}{2}(1 - \alpha)\lambda \sum_{j = 1}^pd_j\beta_j^2\)</span> added to it and the original penalty multiplied by <span class="math inline">\(\alpha\)</span>. For example, <code>grp.mcp.net</code> is the penalty</p>
<p><span class="math display">\[\lambda \sum_{k = 1}^G\alpha p^{MCP}_{\lambda d_k,\gamma}(||\boldsymbol\beta_{g_k}||_2) + \frac{1}{2}(1 - \alpha)\lambda \sum_{j = 1}^pd_j\beta_j^2.
\]</span></p>
</div>
<div id="available-model-families" class="section level2">
<h2><span class="header-section-number">4.3</span> Available Model Families</h2>
<p>The following models are available currently.</p>
<table>
<colgroup>
<col width="31%"></col>
<col width="29%"></col>
<col width="39%"></col>
</colgroup>
<thead>
<tr class="header">
<th>Model</th>
<th>Option Name</th>
<th>Loss Form</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td>Linear Regression</td>
<td><code>gaussian</code></td>
<td><span class="math inline">\(\frac{1}{2n}\sum_{i=1}^n(y_i - x_i^T\beta) ^ 2\)</span></td>
</tr>
<tr class="even">
<td>Logistic Regression</td>
<td><code>binomial</code></td>
<td><span class="math inline">\(-\frac{1}{n}\sum_{i=1}^n\left[y_i x_i^T\beta - \log (1 + \exp\{ x_i^T\beta \} ) \right]\)</span></td>
</tr>
</tbody>
</table>
<p>There are plans to include support for multiple responses, binomial models (not just logistic regression), Cox’s proportional hazards model, and more if requested.</p>
</div>
</div>
<div id="fitting-multiple-penalties-at-once" class="section level1">
<h1><span class="header-section-number">5</span> Fitting multiple penalties at once</h1>
<p>The oem algorithm is well-suited to quickly estimate a solution path for multiple penalties simultaneously if the number of variables is not too large. The oem algorithm is only efficient for multiple penalties for linear models.</p>
<p>For the group lasso penalty, the <code>groups</code> argument must be used. <code>groups</code> should be a vector which indicates the group number for each variable.</p>
<div class="sourceCode" id="cb7"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb7-1" data-line-number="1">fit2 &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> x, <span class="dt">y =</span> y, <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;mcp&quot;</span>, <span class="st">&quot;grp.lasso&quot;</span>, <span class="st">&quot;grp.mcp&quot;</span>),</a>
<a class="sourceLine" id="cb7-2" data-line-number="2">            <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>))</a></code></pre></div>
<p>Plot the solution paths for all models</p>
<div class="sourceCode" id="cb8"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb8-1" data-line-number="1"><span class="kw">layout</span>(<span class="kw">matrix</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">4</span>, <span class="dt">ncol =</span> <span class="dv">2</span>))</a>
<a class="sourceLine" id="cb8-2" data-line-number="2"><span class="kw">plot</span>(fit2, <span class="dt">which.model =</span> <span class="dv">1</span>, <span class="dt">xvar =</span> <span class="st">&quot;lambda&quot;</span>)</a>
<a class="sourceLine" id="cb8-3" data-line-number="3"><span class="kw">plot</span>(fit2, <span class="dt">which.model =</span> <span class="dv">2</span>, <span class="dt">xvar =</span> <span class="st">&quot;lambda&quot;</span>)</a>
<a class="sourceLine" id="cb8-4" data-line-number="4"><span class="kw">plot</span>(fit2, <span class="dt">which.model =</span> <span class="dv">3</span>, <span class="dt">xvar =</span> <span class="st">&quot;lambda&quot;</span>)</a>
<a class="sourceLine" id="cb8-5" data-line-number="5"><span class="kw">plot</span>(fit2, <span class="dt">which.model =</span> <span class="st">&quot;grp.mcp&quot;</span>, <span class="dt">xvar =</span> <span class="st">&quot;lambda&quot;</span>)</a></code></pre></div>
<p><img src="" /><!-- --></p>
<div id="timing-comparison" class="section level2">
<h2><span class="header-section-number">5.1</span> Timing Comparison</h2>
<p>The following is a demonstration of oem’s efficiency for computing solutions for tuning parameter paths for multiple penalties at once.</p>
<div id="linear-regression" class="section level3">
<h3><span class="header-section-number">5.1.1</span> Linear Regression</h3>
<p>The efficiency oem for fitting multiple penalties at once only applies to linear models. However, for linear models it is quite efficient, even for a high number of tuning parameters for many different penalties.</p>
<div class="sourceCode" id="cb9"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb9-1" data-line-number="1">nobs  &lt;-<span class="st"> </span><span class="fl">1e5</span></a>
<a class="sourceLine" id="cb9-2" data-line-number="2">nvars &lt;-<span class="st"> </span><span class="dv">100</span></a>
<a class="sourceLine" id="cb9-3" data-line-number="3">x2 &lt;-<span class="st"> </span><span class="kw">matrix</span>(<span class="kw">rnorm</span>(nobs <span class="op">*</span><span class="st"> </span>nvars), <span class="dt">ncol =</span> nvars)</a>
<a class="sourceLine" id="cb9-4" data-line-number="4">y2 &lt;-<span class="st"> </span><span class="kw">drop</span>(x2 <span class="op">%*%</span><span class="st"> </span><span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">-0.5</span>, <span class="fl">-0.5</span>, <span class="dv">1</span>, <span class="kw">rep</span>(<span class="dv">0</span>, nvars <span class="op">-</span><span class="st"> </span><span class="dv">5</span>))) <span class="op">+</span><span class="st"> </span><span class="kw">rnorm</span>(nobs, <span class="dt">sd =</span> <span class="dv">4</span>)</a>
<a class="sourceLine" id="cb9-5" data-line-number="5"></a>
<a class="sourceLine" id="cb9-6" data-line-number="6"><span class="kw">system.time</span>(fit2a &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> x2, <span class="dt">y =</span> y2, <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;grp.lasso&quot;</span>),</a>
<a class="sourceLine" id="cb9-7" data-line-number="7">                         <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>), <span class="dt">nlambda =</span> 100L))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    0.23    0.00    0.23</code></pre>
<div class="sourceCode" id="cb11"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb11-1" data-line-number="1"><span class="kw">system.time</span>(fit2b &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> x2, <span class="dt">y =</span> y2, </a>
<a class="sourceLine" id="cb11-2" data-line-number="2">                         <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;grp.lasso&quot;</span>, <span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;mcp&quot;</span>, </a>
<a class="sourceLine" id="cb11-3" data-line-number="3">                                     <span class="st">&quot;scad&quot;</span>, <span class="st">&quot;elastic.net&quot;</span>, <span class="st">&quot;grp.mcp&quot;</span>,</a>
<a class="sourceLine" id="cb11-4" data-line-number="4">                                     <span class="st">&quot;grp.scad&quot;</span>, <span class="st">&quot;sparse.grp.lasso&quot;</span>),</a>
<a class="sourceLine" id="cb11-5" data-line-number="5">                         <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>), <span class="dt">nlambda =</span> 100L))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    0.22    0.03    0.25</code></pre>
<div class="sourceCode" id="cb13"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb13-1" data-line-number="1"><span class="kw">system.time</span>(fit2c &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> x2, <span class="dt">y =</span> y2, </a>
<a class="sourceLine" id="cb13-2" data-line-number="2">                         <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;grp.lasso&quot;</span>, <span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;mcp&quot;</span>, </a>
<a class="sourceLine" id="cb13-3" data-line-number="3">                                     <span class="st">&quot;scad&quot;</span>, <span class="st">&quot;elastic.net&quot;</span>, <span class="st">&quot;grp.mcp&quot;</span>,</a>
<a class="sourceLine" id="cb13-4" data-line-number="4">                                     <span class="st">&quot;grp.scad&quot;</span>, <span class="st">&quot;sparse.grp.lasso&quot;</span>),</a>
<a class="sourceLine" id="cb13-5" data-line-number="5">                         <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>), <span class="dt">nlambda =</span> 500L))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    0.34    0.01    0.39</code></pre>
</div>
<div id="logistic-regression" class="section level3">
<h3><span class="header-section-number">5.1.2</span> Logistic Regression</h3>
<p>It is still more efficient to fit multiple penalties at once instead of individually for logistic regression, but the benefit is not as dramatic as for linear models.</p>
<div class="sourceCode" id="cb15"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb15-1" data-line-number="1">nobs  &lt;-<span class="st"> </span><span class="fl">5e4</span></a>
<a class="sourceLine" id="cb15-2" data-line-number="2">nvars &lt;-<span class="st"> </span><span class="dv">100</span></a>
<a class="sourceLine" id="cb15-3" data-line-number="3">x2 &lt;-<span class="st"> </span><span class="kw">matrix</span>(<span class="kw">rnorm</span>(nobs <span class="op">*</span><span class="st"> </span>nvars), <span class="dt">ncol =</span> nvars)</a>
<a class="sourceLine" id="cb15-4" data-line-number="4"></a>
<a class="sourceLine" id="cb15-5" data-line-number="5">y2 &lt;-<span class="st"> </span><span class="kw">rbinom</span>(nobs, <span class="dv">1</span>, <span class="dt">prob =</span> <span class="dv">1</span> <span class="op">/</span><span class="st"> </span>(<span class="dv">1</span> <span class="op">+</span><span class="st"> </span><span class="kw">exp</span>(<span class="op">-</span><span class="kw">drop</span>(x2 <span class="op">%*%</span><span class="st"> </span><span class="kw">c</span>(<span class="fl">0.15</span>, <span class="fl">0.15</span>, <span class="fl">-0.15</span>, <span class="fl">-0.15</span>, <span class="fl">0.25</span>, <span class="kw">rep</span>(<span class="dv">0</span>, nvars <span class="op">-</span><span class="st"> </span><span class="dv">5</span>))))))</a>
<a class="sourceLine" id="cb15-6" data-line-number="6"></a>
<a class="sourceLine" id="cb15-7" data-line-number="7"></a>
<a class="sourceLine" id="cb15-8" data-line-number="8"><span class="kw">system.time</span>(fit2a &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> x2, <span class="dt">y =</span> y2, <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;grp.lasso&quot;</span>),</a>
<a class="sourceLine" id="cb15-9" data-line-number="9">                         <span class="dt">family =</span> <span class="st">&quot;binomial&quot;</span>,</a>
<a class="sourceLine" id="cb15-10" data-line-number="10">                         <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>), <span class="dt">nlambda =</span> 100L))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    2.59    0.02    2.64</code></pre>
<div class="sourceCode" id="cb17"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb17-1" data-line-number="1"><span class="kw">system.time</span>(fit2b &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> x2, <span class="dt">y =</span> y2, <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;grp.lasso&quot;</span>, <span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;mcp&quot;</span>, <span class="st">&quot;scad&quot;</span>, <span class="st">&quot;elastic.net&quot;</span>),</a>
<a class="sourceLine" id="cb17-2" data-line-number="2">                         <span class="dt">family =</span> <span class="st">&quot;binomial&quot;</span>,</a>
<a class="sourceLine" id="cb17-3" data-line-number="3">                         <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>), <span class="dt">nlambda =</span> 100L))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##   10.53    0.14   10.83</code></pre>
</div>
</div>
</div>
<div id="cross-validation" class="section level1">
<h1><span class="header-section-number">6</span> Cross Validation</h1>
<p>Here we use the <code>nfolds</code> argument to specify the number of folds for <span class="math inline">\(k\)</span>-fold cross validation</p>
<div class="sourceCode" id="cb19"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb19-1" data-line-number="1"><span class="kw">system.time</span>(cvfit1 &lt;-<span class="st"> </span><span class="kw">cv.oem</span>(<span class="dt">x =</span> x, <span class="dt">y =</span> y, </a>
<a class="sourceLine" id="cb19-2" data-line-number="2">                             <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;mcp&quot;</span>, </a>
<a class="sourceLine" id="cb19-3" data-line-number="3">                                         <span class="st">&quot;grp.lasso&quot;</span>, <span class="st">&quot;grp.mcp&quot;</span>), </a>
<a class="sourceLine" id="cb19-4" data-line-number="4">                             <span class="dt">gamma =</span> <span class="dv">2</span>,</a>
<a class="sourceLine" id="cb19-5" data-line-number="5">                             <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>), </a>
<a class="sourceLine" id="cb19-6" data-line-number="6">                             <span class="dt">nfolds =</span> <span class="dv">10</span>))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    1.53    0.14    1.67</code></pre>
<p>Plot the cross validation mean squared error results for each model</p>
<div class="sourceCode" id="cb21"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb21-1" data-line-number="1"><span class="kw">layout</span>(<span class="kw">matrix</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">4</span>, <span class="dt">ncol =</span> <span class="dv">2</span>))</a>
<a class="sourceLine" id="cb21-2" data-line-number="2"><span class="kw">plot</span>(cvfit1, <span class="dt">which.model =</span> <span class="dv">1</span>)</a>
<a class="sourceLine" id="cb21-3" data-line-number="3"><span class="kw">plot</span>(cvfit1, <span class="dt">which.model =</span> <span class="dv">2</span>)</a>
<a class="sourceLine" id="cb21-4" data-line-number="4"><span class="kw">plot</span>(cvfit1, <span class="dt">which.model =</span> <span class="dv">3</span>)</a>
<a class="sourceLine" id="cb21-5" data-line-number="5"><span class="kw">plot</span>(cvfit1, <span class="dt">which.model =</span> <span class="dv">4</span>)</a></code></pre></div>
<p><img src="" /><!-- --></p>
<div id="extremely-fast-cross-validation-for-linear-models" class="section level2">
<h2><span class="header-section-number">6.1</span> Extremely Fast Cross Validation for Linear Models</h2>
<p>The function <code>xval.oem</code> offers accelerated cross validation for penalized linear models. In many cases is is orders of magnitude faster than cv.oem. It is only recommended for scenarios where the number of observations is larger than the number of variables. In addition to the computational gains in single-core usage, it also benefits from parallelizaton using OpenMP (instead of using foreach, as used by cv.oem). For large enough problems, it has on a similar order of computation time as just fitting one OEM model.</p>
<div class="sourceCode" id="cb22"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb22-1" data-line-number="1">nobsc  &lt;-<span class="st"> </span><span class="fl">1e5</span></a>
<a class="sourceLine" id="cb22-2" data-line-number="2">nvarsc &lt;-<span class="st"> </span><span class="dv">100</span></a>
<a class="sourceLine" id="cb22-3" data-line-number="3">xc &lt;-<span class="st"> </span><span class="kw">matrix</span>(<span class="kw">rnorm</span>(nobsc <span class="op">*</span><span class="st"> </span>nvarsc), <span class="dt">ncol =</span> nvarsc)</a>
<a class="sourceLine" id="cb22-4" data-line-number="4">yc &lt;-<span class="st"> </span><span class="kw">drop</span>(xc <span class="op">%*%</span><span class="st"> </span><span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">-0.5</span>, <span class="fl">-0.5</span>, <span class="dv">1</span>, <span class="kw">rep</span>(<span class="dv">0</span>, nvarsc <span class="op">-</span><span class="st"> </span><span class="dv">5</span>))) <span class="op">+</span><span class="st"> </span><span class="kw">rnorm</span>(nobsc, <span class="dt">sd =</span> <span class="dv">4</span>)</a>
<a class="sourceLine" id="cb22-5" data-line-number="5"></a>
<a class="sourceLine" id="cb22-6" data-line-number="6"><span class="kw">system.time</span>(cvalfit1 &lt;-<span class="st"> </span><span class="kw">cv.oem</span>(<span class="dt">x =</span> xc, <span class="dt">y =</span> yc, <span class="dt">penalty =</span> <span class="st">&quot;lasso&quot;</span>, </a>
<a class="sourceLine" id="cb22-7" data-line-number="7">                               <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>), </a>
<a class="sourceLine" id="cb22-8" data-line-number="8">                               <span class="dt">nfolds =</span> <span class="dv">10</span>))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    5.30    0.89    6.27</code></pre>
<div class="sourceCode" id="cb24"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb24-1" data-line-number="1"><span class="kw">system.time</span>(xvalfit1 &lt;-<span class="st"> </span><span class="kw">xval.oem</span>(<span class="dt">x =</span> xc, <span class="dt">y =</span> yc, <span class="dt">penalty =</span> <span class="st">&quot;lasso&quot;</span>,</a>
<a class="sourceLine" id="cb24-2" data-line-number="2">                                 <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>), </a>
<a class="sourceLine" id="cb24-3" data-line-number="3">                                 <span class="dt">nfolds =</span> <span class="dv">10</span>))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    0.86    0.08    0.94</code></pre>
<div class="sourceCode" id="cb26"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb26-1" data-line-number="1"><span class="kw">system.time</span>(xvalfit2 &lt;-<span class="st"> </span><span class="kw">xval.oem</span>(<span class="dt">x =</span> xc, <span class="dt">y =</span> yc, <span class="dt">penalty =</span> <span class="st">&quot;lasso&quot;</span>,</a>
<a class="sourceLine" id="cb26-2" data-line-number="2">                                 <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>), </a>
<a class="sourceLine" id="cb26-3" data-line-number="3">                                 <span class="dt">nfolds =</span> <span class="dv">10</span>, <span class="dt">ncores =</span> <span class="dv">2</span>))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    1.03    0.03    0.69</code></pre>
<div class="sourceCode" id="cb28"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb28-1" data-line-number="1"><span class="kw">system.time</span>(ofit1 &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> xc, <span class="dt">y =</span> yc, <span class="dt">penalty =</span> <span class="st">&quot;lasso&quot;</span>,</a>
<a class="sourceLine" id="cb28-2" data-line-number="2">                         <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>)))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    0.20    0.02    0.22</code></pre>
</div>
<div id="evaluation-metrics" class="section level2">
<h2><span class="header-section-number">6.2</span> Evaluation Metrics</h2>
<p>A variety of evaluation metrics can be used for cross validation. The available metrics can be found in the table below</p>
<table>
<thead>
<tr class="header">
<th>Model</th>
<th>Metric</th>
<th><code>type.measure=</code></th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td><strong>Linear Regression</strong></td>
<td>Mean squared error</td>
<td><code>mse</code> or <code>deviance</code></td>
</tr>
<tr class="even">
<td></td>
<td>Mean absolute error</td>
<td><code>mae</code></td>
</tr>
<tr class="odd">
<td>———————</td>
<td>———————————-</td>
<td>————————</td>
</tr>
<tr class="even">
<td><strong>Logistic Regression</strong></td>
<td>Deviance</td>
<td><code>deviance</code></td>
</tr>
<tr class="odd">
<td></td>
<td>Area under the ROC curve</td>
<td><code>auc</code></td>
</tr>
<tr class="even">
<td></td>
<td>Misclassification Rate</td>
<td><code>class</code></td>
</tr>
<tr class="odd">
<td></td>
<td>Mean squared error of fitted mean</td>
<td><code>mse</code></td>
</tr>
<tr class="even">
<td></td>
<td>Mean absolute error of fitted mean</td>
<td><code>mae</code></td>
</tr>
</tbody>
</table>
<p>Consider a binary outcome setting with logistic regression.</p>
<div class="sourceCode" id="cb30"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb30-1" data-line-number="1">nobs  &lt;-<span class="st"> </span><span class="fl">2e3</span></a>
<a class="sourceLine" id="cb30-2" data-line-number="2">nvars &lt;-<span class="st"> </span><span class="dv">20</span></a>
<a class="sourceLine" id="cb30-3" data-line-number="3">x &lt;-<span class="st"> </span><span class="kw">matrix</span>(<span class="kw">runif</span>(nobs <span class="op">*</span><span class="st"> </span>nvars, <span class="dt">max =</span> <span class="dv">2</span>), <span class="dt">ncol =</span> nvars)</a>
<a class="sourceLine" id="cb30-4" data-line-number="4"></a>
<a class="sourceLine" id="cb30-5" data-line-number="5">y &lt;-<span class="st"> </span><span class="kw">rbinom</span>(nobs, <span class="dv">1</span>, <span class="dt">prob =</span> <span class="dv">1</span> <span class="op">/</span><span class="st"> </span>(<span class="dv">1</span> <span class="op">+</span><span class="st"> </span><span class="kw">exp</span>(<span class="op">-</span><span class="kw">drop</span>(x <span class="op">%*%</span><span class="st"> </span><span class="kw">c</span>(<span class="fl">0.25</span>, <span class="dv">-1</span>, <span class="dv">-1</span>, <span class="fl">-0.5</span>, <span class="fl">-0.5</span>, <span class="fl">-0.25</span>, <span class="kw">rep</span>(<span class="dv">0</span>, nvars <span class="op">-</span><span class="st"> </span><span class="dv">6</span>))))))</a></code></pre></div>
<div id="misclassification-rate" class="section level3">
<h3><span class="header-section-number">6.2.1</span> Misclassification Rate</h3>
<div class="sourceCode" id="cb31"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb31-1" data-line-number="1">cvfit2 &lt;-<span class="st"> </span><span class="kw">cv.oem</span>(<span class="dt">x =</span> x, <span class="dt">y =</span> y, <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;mcp&quot;</span>, </a>
<a class="sourceLine" id="cb31-2" data-line-number="2">                                           <span class="st">&quot;grp.lasso&quot;</span>, <span class="st">&quot;grp.mcp&quot;</span>), </a>
<a class="sourceLine" id="cb31-3" data-line-number="3">                 <span class="dt">family =</span> <span class="st">&quot;binomial&quot;</span>,</a>
<a class="sourceLine" id="cb31-4" data-line-number="4">                 <span class="dt">type.measure =</span> <span class="st">&quot;class&quot;</span>,</a>
<a class="sourceLine" id="cb31-5" data-line-number="5">                 <span class="dt">gamma =</span> <span class="dv">2</span>,</a>
<a class="sourceLine" id="cb31-6" data-line-number="6">                 <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">10</span>, <span class="dt">each =</span> <span class="dv">2</span>), </a>
<a class="sourceLine" id="cb31-7" data-line-number="7">                 <span class="dt">nfolds =</span> <span class="dv">10</span>, <span class="dt">standardize =</span> <span class="ot">FALSE</span>)</a></code></pre></div>
<p><img src="" /><!-- --></p>
<p>In this case, misclassification rate is not the best indicator of performance. The classes here are imbalanced:</p>
<div class="sourceCode" id="cb32"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb32-1" data-line-number="1"><span class="kw">mean</span>(y)</a></code></pre></div>
<pre><code>## [1] 0.0705</code></pre>
</div>
<div id="area-under-the-roc-curve" class="section level3">
<h3><span class="header-section-number">6.2.2</span> Area Under the ROC Curve</h3>
<p>Area under the ROC curve is an alternative classification metric to misclassification rate. It is available by setting <code>type.measure = &quot;auc&quot;</code>.</p>
<div class="sourceCode" id="cb34"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb34-1" data-line-number="1">cvfit2 &lt;-<span class="st"> </span><span class="kw">cv.oem</span>(<span class="dt">x =</span> x, <span class="dt">y =</span> y, <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;mcp&quot;</span>, </a>
<a class="sourceLine" id="cb34-2" data-line-number="2">                                           <span class="st">&quot;grp.lasso&quot;</span>, <span class="st">&quot;grp.mcp&quot;</span>), </a>
<a class="sourceLine" id="cb34-3" data-line-number="3">                 <span class="dt">family =</span> <span class="st">&quot;binomial&quot;</span>,</a>
<a class="sourceLine" id="cb34-4" data-line-number="4">                 <span class="dt">type.measure =</span> <span class="st">&quot;auc&quot;</span>,</a>
<a class="sourceLine" id="cb34-5" data-line-number="5">                 <span class="dt">gamma =</span> <span class="dv">2</span>,</a>
<a class="sourceLine" id="cb34-6" data-line-number="6">                 <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">10</span>, <span class="dt">each =</span> <span class="dv">2</span>), </a>
<a class="sourceLine" id="cb34-7" data-line-number="7">                 <span class="dt">nfolds =</span> <span class="dv">10</span>, <span class="dt">standardize =</span> <span class="ot">FALSE</span>)</a></code></pre></div>
<p><img src="" /><!-- --></p>
</div>
</div>
</div>
<div id="methods-for-very-large-scale-problems" class="section level1">
<h1><span class="header-section-number">7</span> Methods for Very Large Scale Problems</h1>
<div id="oem-with-precomputed-xtx-xty-for-linear-models" class="section level2">
<h2><span class="header-section-number">7.1</span> OEM with Precomputed <span class="math inline">\(X^TX\)</span>, <span class="math inline">\(X^TY\)</span> for Linear Models</h2>
<p>With a very large dataset and computing cluster, the total size of a dataset may be very large, but if the number of variables is only moderately large (on the order of a few thousands) <span class="math inline">\(X^TX\)</span> and <span class="math inline">\(X^TY\)</span> may not be large and may already be available from other computations or can be computed trivially in parallel. The function <code>oem.xtx</code> computes penalized linear regression models using the OEM algorithm only using <span class="math inline">\(X^TX\)</span> and <span class="math inline">\(X^TY\)</span>. Standardization can be achieved by providing a vector of scaling factors (usually the standard deviations of the columns of x). The function is used like the following:</p>
<div class="sourceCode" id="cb35"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb35-1" data-line-number="1">xtx &lt;-<span class="st"> </span><span class="kw">crossprod</span>(xc) <span class="op">/</span><span class="st"> </span><span class="kw">nrow</span>(xc)</a>
<a class="sourceLine" id="cb35-2" data-line-number="2">xty &lt;-<span class="st"> </span><span class="kw">crossprod</span>(xc, yc) <span class="op">/</span><span class="st"> </span><span class="kw">nrow</span>(xc)</a>
<a class="sourceLine" id="cb35-3" data-line-number="3"></a>
<a class="sourceLine" id="cb35-4" data-line-number="4"></a>
<a class="sourceLine" id="cb35-5" data-line-number="5"><span class="kw">system.time</span>(fit &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> xc, <span class="dt">y =</span> yc, </a>
<a class="sourceLine" id="cb35-6" data-line-number="6">                       <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;grp.lasso&quot;</span>), </a>
<a class="sourceLine" id="cb35-7" data-line-number="7">                       <span class="dt">standardize =</span> <span class="ot">FALSE</span>, <span class="dt">intercept =</span> <span class="ot">FALSE</span>,</a>
<a class="sourceLine" id="cb35-8" data-line-number="8">                       <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>)))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    0.19    0.04    0.23</code></pre>
<div class="sourceCode" id="cb37"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb37-1" data-line-number="1"><span class="kw">system.time</span>(fit.xtx &lt;-<span class="st"> </span><span class="kw">oem.xtx</span>(<span class="dt">xtx =</span> xtx, <span class="dt">xty =</span> xty, </a>
<a class="sourceLine" id="cb37-2" data-line-number="2">                               <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;grp.lasso&quot;</span>), </a>
<a class="sourceLine" id="cb37-3" data-line-number="3">                               <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>))  )  </a></code></pre></div>
<pre><code>##    user  system elapsed 
##    0.01    0.00    0.02</code></pre>
<div class="sourceCode" id="cb39"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb39-1" data-line-number="1"><span class="kw">max</span>(<span class="kw">abs</span>(fit<span class="op">$</span>beta[[<span class="dv">1</span>]][<span class="op">-</span><span class="dv">1</span>,] <span class="op">-</span><span class="st"> </span>fit.xtx<span class="op">$</span>beta[[<span class="dv">1</span>]]))</a></code></pre></div>
<pre><code>## [1] 9.992007e-15</code></pre>
<div class="sourceCode" id="cb41"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb41-1" data-line-number="1"><span class="kw">max</span>(<span class="kw">abs</span>(fit<span class="op">$</span>beta[[<span class="dv">2</span>]][<span class="op">-</span><span class="dv">1</span>,] <span class="op">-</span><span class="st"> </span>fit.xtx<span class="op">$</span>beta[[<span class="dv">2</span>]])) </a></code></pre></div>
<pre><code>## [1] 9.992007e-15</code></pre>
<div class="sourceCode" id="cb43"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb43-1" data-line-number="1">col.std &lt;-<span class="st"> </span><span class="kw">apply</span>(xc, <span class="dv">2</span>, sd)</a>
<a class="sourceLine" id="cb43-2" data-line-number="2">fit.xtx.s &lt;-<span class="st"> </span><span class="kw">oem.xtx</span>(<span class="dt">xtx =</span> xtx, <span class="dt">xty =</span> xty, </a>
<a class="sourceLine" id="cb43-3" data-line-number="3">                     <span class="dt">scale.factor =</span> col.std,</a>
<a class="sourceLine" id="cb43-4" data-line-number="4">                     <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;grp.lasso&quot;</span>), </a>
<a class="sourceLine" id="cb43-5" data-line-number="5">                     <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>))  </a></code></pre></div>
</div>
<div id="out-of-memory-computation" class="section level2">
<h2><span class="header-section-number">7.2</span> Out-of-memory Computation</h2>
<p>The OEM package also provides functionality for on-disk computation with the <code>big.oem</code> function, allowing for fitting penalized regression models on datasets too large to fit in memory. The <code>big.oem</code> function uses the tools provided by the <code>bigmemory</code> package, so a big.matrix object must be used for the design matrix.</p>
<div class="sourceCode" id="cb44"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb44-1" data-line-number="1"><span class="kw">set.seed</span>(<span class="dv">123</span>)</a>
<a class="sourceLine" id="cb44-2" data-line-number="2">nrows &lt;-<span class="st"> </span><span class="dv">50000</span></a>
<a class="sourceLine" id="cb44-3" data-line-number="3">ncols &lt;-<span class="st"> </span><span class="dv">100</span></a>
<a class="sourceLine" id="cb44-4" data-line-number="4">bkFile &lt;-<span class="st"> &quot;bigmat.bk&quot;</span></a>
<a class="sourceLine" id="cb44-5" data-line-number="5">descFile &lt;-<span class="st"> &quot;bigmatk.desc&quot;</span></a>
<a class="sourceLine" id="cb44-6" data-line-number="6">bigmat &lt;-<span class="st"> </span><span class="kw">filebacked.big.matrix</span>(<span class="dt">nrow=</span>nrows, <span class="dt">ncol=</span>ncols, <span class="dt">type=</span><span class="st">&quot;double&quot;</span>,</a>
<a class="sourceLine" id="cb44-7" data-line-number="7">                                <span class="dt">backingfile=</span>bkFile, <span class="dt">backingpath=</span><span class="st">&quot;.&quot;</span>,</a>
<a class="sourceLine" id="cb44-8" data-line-number="8">                                <span class="dt">descriptorfile=</span>descFile,</a>
<a class="sourceLine" id="cb44-9" data-line-number="9">                                <span class="dt">dimnames=</span><span class="kw">c</span>(<span class="ot">NULL</span>,<span class="ot">NULL</span>))</a>
<a class="sourceLine" id="cb44-10" data-line-number="10"></a>
<a class="sourceLine" id="cb44-11" data-line-number="11"><span class="co"># Each column value with be the column number multiplied by</span></a>
<a class="sourceLine" id="cb44-12" data-line-number="12"><span class="co"># samples from a standard normal distribution.</span></a>
<a class="sourceLine" id="cb44-13" data-line-number="13"><span class="kw">set.seed</span>(<span class="dv">123</span>)</a>
<a class="sourceLine" id="cb44-14" data-line-number="14"><span class="cf">for</span> (i <span class="cf">in</span> <span class="dv">1</span><span class="op">:</span>ncols) bigmat[,i] =<span class="st"> </span><span class="kw">rnorm</span>(nrows)<span class="op">*</span>i</a>
<a class="sourceLine" id="cb44-15" data-line-number="15"></a>
<a class="sourceLine" id="cb44-16" data-line-number="16">yb &lt;-<span class="st"> </span><span class="kw">rnorm</span>(nrows) <span class="op">+</span><span class="st"> </span>bigmat[,<span class="dv">1</span>] <span class="op">-</span><span class="st"> </span>bigmat[,<span class="dv">2</span>]</a>
<a class="sourceLine" id="cb44-17" data-line-number="17"></a>
<a class="sourceLine" id="cb44-18" data-line-number="18">## out-of-memory computation</a>
<a class="sourceLine" id="cb44-19" data-line-number="19">fit &lt;-<span class="st"> </span><span class="kw">big.oem</span>(<span class="dt">x =</span> bigmat, <span class="dt">y =</span> yb, </a>
<a class="sourceLine" id="cb44-20" data-line-number="20">               <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;grp.lasso&quot;</span>), </a>
<a class="sourceLine" id="cb44-21" data-line-number="21">               <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>))</a>
<a class="sourceLine" id="cb44-22" data-line-number="22"></a>
<a class="sourceLine" id="cb44-23" data-line-number="23">## fitting with in-memory computation</a>
<a class="sourceLine" id="cb44-24" data-line-number="24">fit2 &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> bigmat[,], <span class="dt">y =</span> yb, </a>
<a class="sourceLine" id="cb44-25" data-line-number="25">            <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;grp.lasso&quot;</span>), </a>
<a class="sourceLine" id="cb44-26" data-line-number="26">            <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">5</span>))   </a>
<a class="sourceLine" id="cb44-27" data-line-number="27">           </a>
<a class="sourceLine" id="cb44-28" data-line-number="28"><span class="kw">max</span>(<span class="kw">abs</span>(fit<span class="op">$</span>beta[[<span class="dv">1</span>]] <span class="op">-</span><span class="st"> </span>fit2<span class="op">$</span>beta[[<span class="dv">1</span>]]))            </a></code></pre></div>
<pre><code>## [1] 1.534783e-05</code></pre>
</div>
</div>
<div id="other-features" class="section level1">
<h1><span class="header-section-number">8</span> Other Features</h1>
<div id="parallelization-via-openmp" class="section level2">
<h2><span class="header-section-number">8.1</span> Parallelization via OpenMP</h2>
<p>Computational time can be reduced a little via OpenMP parallelization of the key computational steps of the OEM algorithm. Simply use the <code>ncores</code> argument to access parallelization. There is no need for the foreach package.</p>
<div class="sourceCode" id="cb46"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb46-1" data-line-number="1">nobsc  &lt;-<span class="st"> </span><span class="fl">1e5</span></a>
<a class="sourceLine" id="cb46-2" data-line-number="2">nvarsc &lt;-<span class="st"> </span><span class="dv">500</span></a>
<a class="sourceLine" id="cb46-3" data-line-number="3">xc &lt;-<span class="st"> </span><span class="kw">matrix</span>(<span class="kw">rnorm</span>(nobsc <span class="op">*</span><span class="st"> </span>nvarsc), <span class="dt">ncol =</span> nvarsc)</a>
<a class="sourceLine" id="cb46-4" data-line-number="4">yc &lt;-<span class="st"> </span><span class="kw">drop</span>(xc <span class="op">%*%</span><span class="st"> </span><span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">-0.5</span>, <span class="fl">-0.5</span>, <span class="dv">1</span>, <span class="kw">rep</span>(<span class="dv">0</span>, nvarsc <span class="op">-</span><span class="st"> </span><span class="dv">5</span>))) <span class="op">+</span><span class="st"> </span><span class="kw">rnorm</span>(nobsc, <span class="dt">sd =</span> <span class="dv">4</span>)</a>
<a class="sourceLine" id="cb46-5" data-line-number="5"></a>
<a class="sourceLine" id="cb46-6" data-line-number="6"></a>
<a class="sourceLine" id="cb46-7" data-line-number="7"><span class="kw">system.time</span>(fit &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> xc, <span class="dt">y =</span> yc, </a>
<a class="sourceLine" id="cb46-8" data-line-number="8">                       <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;grp.lasso&quot;</span>), </a>
<a class="sourceLine" id="cb46-9" data-line-number="9">                       <span class="dt">standardize =</span> <span class="ot">FALSE</span>, <span class="dt">intercept =</span> <span class="ot">FALSE</span>,</a>
<a class="sourceLine" id="cb46-10" data-line-number="10">                       <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">25</span>)))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    5.03    0.27    6.12</code></pre>
<div class="sourceCode" id="cb48"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb48-1" data-line-number="1"><span class="kw">system.time</span>(fitp &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> xc, <span class="dt">y =</span> yc, </a>
<a class="sourceLine" id="cb48-2" data-line-number="2">                        <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>, <span class="st">&quot;grp.lasso&quot;</span>), </a>
<a class="sourceLine" id="cb48-3" data-line-number="3">                        <span class="dt">standardize =</span> <span class="ot">FALSE</span>, <span class="dt">intercept =</span> <span class="ot">FALSE</span>,</a>
<a class="sourceLine" id="cb48-4" data-line-number="4">                        <span class="dt">groups =</span> <span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">20</span>, <span class="dt">each =</span> <span class="dv">25</span>), <span class="dt">ncores =</span> <span class="dv">2</span>))</a></code></pre></div>
<pre><code>##    user  system elapsed 
##    5.55    0.12    4.48</code></pre>
</div>
<div id="penalty-adjustment" class="section level2">
<h2><span class="header-section-number">8.2</span> Penalty Adjustment</h2>
<p>If some variables should not be penalized, this can be specified through the use of the <code>penalty.factor</code> argument for all penalties other than the group lasso. For the group lasso, the group-specific weights can be modified by the <code>group.weights</code> argument. <code>penalty.factor</code> should be a vector of length equal to the number of columns in <code>x</code>. Each element in <code>penalty.factor</code> will be multiplied to the applied tuning parameter for each corresponding variable. For example, for a problem with 5 variables (<code>ncol(x) = 5</code>), setting <code>penalty.factor = c(1, 1, 1, 0, 0)</code> will effectively only allow penalization for the first three variables. The <code>group.weights</code> argument should be a vector with length equal to the number of groups. Similarly to <code>penalty.factor</code>, these weights will be multiplied to the penalty applied to each group. <code>penalty.factor</code> and <code>group.weights</code> can also be used to fit the adaptive lasso and adaptive group lasso, respectively.</p>
<p>The following example shows how to fit an adaptive lasso using <code>oem</code></p>
<div class="sourceCode" id="cb50"><pre class="sourceCode r"><code class="sourceCode r"><a class="sourceLine" id="cb50-1" data-line-number="1">nobs  &lt;-<span class="st"> </span><span class="fl">1e4</span></a>
<a class="sourceLine" id="cb50-2" data-line-number="2">nvars &lt;-<span class="st"> </span><span class="dv">102</span></a>
<a class="sourceLine" id="cb50-3" data-line-number="3">x &lt;-<span class="st"> </span><span class="kw">matrix</span>(<span class="kw">rnorm</span>(nobs <span class="op">*</span><span class="st"> </span>nvars), <span class="dt">ncol =</span> nvars)</a>
<a class="sourceLine" id="cb50-4" data-line-number="4">y &lt;-<span class="st"> </span><span class="kw">drop</span>(x <span class="op">%*%</span><span class="st"> </span><span class="kw">c</span>(<span class="fl">0.5</span>, <span class="fl">0.5</span>, <span class="fl">-0.5</span>, <span class="fl">-0.5</span>, <span class="dv">1</span>, <span class="fl">0.5</span>, <span class="kw">rep</span>(<span class="dv">0</span>, nvars <span class="op">-</span><span class="st"> </span><span class="dv">6</span>))) <span class="op">+</span><span class="st"> </span><span class="kw">rnorm</span>(nobs, <span class="dt">sd =</span> <span class="dv">4</span>)</a>
<a class="sourceLine" id="cb50-5" data-line-number="5"></a>
<a class="sourceLine" id="cb50-6" data-line-number="6">lams &lt;-<span class="st"> </span><span class="kw">exp</span>(<span class="kw">seq</span>(<span class="kw">log</span>(<span class="fl">2.5</span>), <span class="kw">log</span>(<span class="fl">5e-5</span>), <span class="dt">length.out =</span> 100L))</a>
<a class="sourceLine" id="cb50-7" data-line-number="7"></a>
<a class="sourceLine" id="cb50-8" data-line-number="8">ols.estimates &lt;-<span class="st"> </span><span class="kw">coef</span>(<span class="kw">lm.fit</span>(<span class="dt">y =</span> y, <span class="dt">x =</span> <span class="kw">cbind</span>(<span class="dv">1</span>, x)))[<span class="op">-</span><span class="dv">1</span>]</a>
<a class="sourceLine" id="cb50-9" data-line-number="9"></a>
<a class="sourceLine" id="cb50-10" data-line-number="10">fit.adaptive &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> x, <span class="dt">y =</span> y, <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;lasso&quot;</span>),</a>
<a class="sourceLine" id="cb50-11" data-line-number="11">                    <span class="dt">penalty.factor =</span> <span class="dv">1</span> <span class="op">/</span><span class="st"> </span><span class="kw">abs</span>(ols.estimates),</a>
<a class="sourceLine" id="cb50-12" data-line-number="12">                    <span class="dt">lambda =</span> lams)</a>
<a class="sourceLine" id="cb50-13" data-line-number="13"></a>
<a class="sourceLine" id="cb50-14" data-line-number="14">group.indicators &lt;-<span class="st"> </span><span class="kw">rep</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">34</span>, <span class="dt">each =</span> <span class="dv">3</span>)</a>
<a class="sourceLine" id="cb50-15" data-line-number="15"></a>
<a class="sourceLine" id="cb50-16" data-line-number="16">## norms of OLS estimates for each group</a>
<a class="sourceLine" id="cb50-17" data-line-number="17">group.norms      &lt;-<span class="st"> </span><span class="kw">sapply</span>(<span class="dv">1</span><span class="op">:</span><span class="dv">34</span>, <span class="cf">function</span>(idx) <span class="kw">sqrt</span>(<span class="kw">sum</span>((ols.estimates[group.indicators <span class="op">==</span><span class="st"> </span>idx]) <span class="op">^</span><span class="st"> </span><span class="dv">2</span>)))</a>
<a class="sourceLine" id="cb50-18" data-line-number="18">fit.adaptive.grp &lt;-<span class="st"> </span><span class="kw">oem</span>(<span class="dt">x =</span> x, <span class="dt">y =</span> y, <span class="dt">penalty =</span> <span class="kw">c</span>(<span class="st">&quot;grp.lasso&quot;</span>),</a>
<a class="sourceLine" id="cb50-19" data-line-number="19">                        <span class="dt">group.weights =</span> <span class="dv">1</span> <span class="op">/</span><span class="st"> </span>group.norms,</a>
<a class="sourceLine" id="cb50-20" data-line-number="20">                        <span class="dt">groups =</span> group.indicators, </a>
<a class="sourceLine" id="cb50-21" data-line-number="21">                        <span class="dt">lambda =</span> lams)</a></code></pre></div>
<p><img src="" /><!-- --></p>
</div>
</div>
<div id="more-information" class="section level1">
<h1><span class="header-section-number">9</span> More Information</h1>
<p>For further information about <code>oem</code>, please visit:</p>
<ul>
<li>The oem site: <a href="http://jaredhuling.github.io/oem">jaredhuling.org/oem</a></li>
<li>The oem source code: <a href="https://github.com/jaredhuling/oem">github.com/jaredhuling/oem</a></li>
<li>The oem paper: <a href="http://www.tandfonline.com/doi/abs/10.1080/00401706.2015.1054436">oem paper</a></li>
</ul>
</div>
<div id="references" class="section level1 unnumbered">
<h1>References</h1>
<div id="refs" class="references">
<div id="ref-xiong16">
<p>Xiong, Shifeng, Bin Dai, Jared Huling, and Peter Z.G. Qian. 2016. “Orthogonalizing EM: A Design-Based Least Squares Algorithm.” <em>Technometrics, in Press</em>.</p>
</div>
</div>
</div>



<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>
